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ABSTRACT 

We introduce a low Mach number equation set for the large-scale numerical 
simulation of carbon-oxygen white dwarfs experiencing a thermonuclear defla- 
gration. Since most of the interesting physics in a Type la supernova transpires 
at Mach numbers from 0.01 to 0.1, such an approach enables both a consider- 
able increase in accuracy and savings in computer time compared with frequently 
used compressible codes. Our equation set is derived from the fully compressible 
equations using low Mach number asymptotics, but without any restriction on 
the size of perturbations in density or temperature. Comparisons with simula- 
tions that use the fully compressible equations validate the low Mach number 
model in regimes where both are applicable. Comparisons to simulations based 
on the more traditional anelastic approximation also demonstrate the agreement 
of these models in the regime for which the anelastic approximation is valid. For 
low Mach number flows with potentially finite amplitude variations in density 
and temperature, the low Mach number model overcomes the limitations of each 
of the more traditional models and can serve as the basis for an accurate and 
efficient simulation tool. 

Subject headings: supernovae: general — white dwarfs — hydrodynamics - 
nuclear reactions, nucleosynthesis, abundances — convection — methods: nu- 
merical 



1. Introduction 

A broad range of interesting phenomena in science and engineering occur in a low 
Mach number regime in which the fluid velocity is much less than the speed of sound. 
Several low Mach number schemes have been developed to exploit this separation of scales; 
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these models capture the fluid dynamics of interest without the need to resolve acoustic 
wave propagation. Physically, one can think of the solution to a low Mach number model 
as supporting infinitely fast acoustic equilibration rather than finite-velocity acoustic wave 
propagation. Mathematically, this is manifest in the addition of a constraint on the velocity 
field to the system of evolution equations. This velocity constraint can be translated into 
an elliptic equation for pressure that expresses the equilibration process. Because explicit 
discretization schemes for the low Mach number system are limited by the fluid velocity 
and not by the sound speed, they often gain several orders of magnitude in computational 
efficiency over the traditional compressible approach. 

The simplest low Mach number model is expressed by the incompressible Navier-Stokes 
equations for a constant density fluid. Generalizations that incorporate variations in den- 
sity include the Boussinesq approximation (Boussinesq 1903), which allows heating- induced 
buoyancy in a constant density background, and the anelastic atmospheric (Batchelor 1953; 
Ogura & Phillips 1962; Dutton & Fichtl 1969; Gough 1969; Lipps & Hemler 1982, 1985; Lipps 
1990; Wilhelmson & Ogura 1972) and stellar (Latour et al. 1976; Gilman & Glatzmaier 1981; 
Glatzmaier 1984) approximations that include the effect of large-scale background stratifica- 
tion in the fluid density and pressure but assume small thermodynamic perturbations from 
the background. Low Mach number models for chemical combustion (Rehm & Baum 1978; 
Majda & Sethian 1985; Day & Bell 2000) and nuclear burning (Bell et al. 2004) incorporate 
large compressibility effects due to chemical/nuclear reactions and thermal processes with a 
spatially constant background pressure. 

Low Mach number models to date have, with one exception, allowed either zero vol- 
umetric changes (incompressible Navier-Stokes, Boussinesq) or changes due only to local 
heating effects (low-speed combustion, nuclear burning), or to large-scale background strat- 
ification (anelastic). The only low Mach number model that incorporates both finite local 
expansion due to heating and volume changes due to background stratification is the pseudo- 
incompressible equation set for the terrestrial atmosphere, introduced by Durran (1989) and 
rigorously derived using low Mach number asymptotics by Botta et al. (1999). The formula- 
tion of the pseudo-incompressible constraint assumes the ideal gas equation of state, which 
results in a simplification of terms that does not hold for more general equations of state 
and non-trivial changes in composition. 

Our anticipated application for this model is the convective, ignition, and early prop- 
agation phases of Type la supernovae, but the model should be applicable to many other 
problems, such as Type I X-ray bursts (see, e.g., recent work by Lin (2005) for an alternate 
form of the low Mach number approach for Type I X-ray bursts), classical novae, and or- 
dinary convection in stars. Events such as Type la supernovae are characterized by a large 
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range in length scales, from the O(10 -4 ) to O(IC) 1 ) cm scale of the flame to the O(10 8 ) cm 
scale of the white dwarf. The range in timescales is equally impressive, from the 100 years 
of convection that precedes ignition to the one second duration of the explosion. Presently, 
most large-scale calculations focus on the explosion itself, beginning with several seeded hot 
spots to begin the runaway. In the last minutes of the convective phase, velocities reach 
~ 1% of the sound speed (Woosley 2001; Woosley et al. 2004), with temperature fluctua- 
tions of at most 5% (Wunsch & Woosley 2004). These speeds are too slow for compressible 
codes to accurately follow. For this reason, simulation of the convective and ignition phases 
has seen only limited numerical work, e.g. Hoflich & Stein (2002). Recent analytic work 
(Woosley et al. 2004) suggests that full star simulations are needed to accurately capture 
the convective flows and yield the spatial, temporal, and size distribution of the hot spots 
that seed the explosion. Three-dimensional anelastic calculations (Kuhlen et al. 2005) have 
shown that a dipole velocity field dominates the convection. 

The low Mach number model presented here, like the anelastic model, will be capable 
of the long time integration necessary to follow the convection. Unlike the anelastic model, 
however, the low Mach number approach continues to be valid as the variation in density and 
temperature increases in the flame bubbles that evolve in the early phases of the explosion. 
Following the evolution from the convection through the early phases of the explosion is the 
eventual target of our low Mach number methodology. 

Although the derivation of the low Mach number equation set will be general, the 
example we will consider is a simplified problem without reactions or thermal conduction, and 
with a time-independent radially symmetric form of self-gravity. The focus of the examples 
is to demonstrate the ability of the low Mach number model to accurately represent the 
hydrodynamics. We will compare the simulations based on the low Mach number approach 
to simulations based on the fully compressible equation set, where applicable, and to the 
traditional anelastic approach, where it is applicable. We will show that the low Mach 
number algorithm works well for very low Mach number flows, with validation presented up 
to Mach numbers ~ 0.2. 

In the next section we derive the low Mach number equations, and in Section 3 we discuss 
the numerical implementation. Section 4 contains numerical comparisons with compressible 
and anelastic simulations, and in the final section we discuss our conclusions and future 
work. 



-4- 



2. Low Mach Number Model 

We begin with the fully compressible equations governing motion in the stellar environ- 
ment as described, for example, in Bell et al. (2004) 

^ + V-( P U) = , (1) 

1 V ■ (pUU) + Vp = -pge r , (2) 



dt 

dpE 

~dT 

dpX k 



V • ( P UE + P U) = V • («VT) - pg(U ■ e r ) - P^k , (3) 

k 

, V ■ (pUX k ) = pu k . (4) 

Here p, U, T, and p are the density, velocity, temperature and pressure, respectively, and 
E — e + U • U/2 is the total energy with e representing the internal energy. In addition, 
X k is the abundance of the kth. isotope, with associated production rate, u k , and energy 
release, q k . Finally, g(r) is the radially dependent gravitational acceleration (resulting from 
spherically symmetric self-gravity), e r is the unit vector in the radial direction, and k is the 
thermal conductivity. The Reynolds number of flows in a typical white dwarf is sufficiently 
large that we neglect viscosity here, though viscous terms could easily be included in the 
model and the numerical methodology. 

For the stellar conditions being considered the pressure contains contributions from ions, 
radiation, and electrons. Thus 

p = p(p, T, X k ) = p ion + p rad + p e i e , (5) 

where 

pksT 4 

Pion = -f , Prad = OT /3 , 

Am p 

and Peie is the contribution to the thermodynamic pressure due to fermions. In these expres- 
sions m p is the mass of the proton, a is related to the Stefan-Boltzmann constant o = ac/4, 
c is the speed of light, A = ^2 k X k A k , A k is the atomic number of the kth. isotope, and ks is 
Boltzmann's constant. The ionic component has the form associated with an ideal gas but 
the radiation and electron pressure components do not. We use a stellar equation of state 
as implemented in Timmes & Swesty (2000). 

As a prelude to developing the low Mach number equations, we first rewrite the energy 
equation (Eq. [3]) in terms of the enthalpy, h = e + pj p, 

9 ~Dt ~~Dt = V ' ^ ^ ~ ^ PqhUJh = pH ' ^ 

k 
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where we introduce H to represent the enthalpy source terms. 

Our goal in this section is to derive a model for low speed flows in a hydrostatically 
balanced, radially stratified, background that removes acoustic waves yet allows for the 
development of finite amplitude temperature and density variations. We thus posit the 
existence of a background state with pressure, density and temperature po(r, t), po(r, t) and 
T (r,t) satisfying both the equation of state and hydrostatic equilibrium. Because we will 
neglect reaction terms that could potentially alter the large-scale pressure distribution within 
the star, for the purposes of this paper we will neglect time variation of the background state, 
i.e., we will assume dpo/dt = dpo/dt = dT /dt = 0. 

In order to understand the behavior of the system, we examine the balance of terms as a 
function of Mach number, M = U/c s (c s is the speed of sound), which will be assumed to be 
small. We re-write the momentum equation in nondimensional coordinates, where the space 
and time coordinates as well as the density, velocity and pressure, are scaled by characteristic 
values L re f, t Te {, p rc f,Prcf, and U Te {, respectively. For the problem scale of interest U ve { is a 
typical advective velocity, L re{ a typical length scale, t Te f = L ref /?7 re f, and p TC f = p re f , 
where c Sic{ is a characteristic value of c s . We define a scaling for g in terms of the pressure 
scale height, H rci = p Tci /(p rci g). 

The momentum equation in nondimensional coordinates (t = t/t ie f, etc.), and exploiting 
hydrostatic equilibrium of the reference state, has the form 

a -§ + v . am + ^v ( p - « = -^jp - , 

For the large-scale near-equilibrium behavior, we set L re f = H re {, getting 

^ + V • (pUU) + ^V(p - Po) = ~{p - p~ )ge r , 

where g = g / (p rc f/(PvcfH rci )). 

Since all nondimensional terms are 0(1), it is clear that to maintain a long-term balance, 
both (p — po) and (p — p~ ) must be 0(M 2 ). This is consistent with the traditional anelastic 
approximation; once the density perturbation is assumed small the approximation V-(pqU) = 
follows from the continuity equation, and a linearized temperature-density relationship 
can be used to replace the buoyancy term in the momentum equation by one dependent on 
temperature (or entropy) rather than density. 

Here, however, we are interested in finite amplitude density perturbations. In this case, 
it is possible for the model to breakdown in long time integrations should the flow accelerate 
to the point that M is no longer small. We would, nevertheless, expect the low Mach number 
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model to remain valid for a limited period of time. The behavior of the model for finite time 
intervals can be examined by considering a shorter time scale, t re f = U rc f/g, defined such 
that on this time scale the buoyancy forcing from finite amplitude density perturbations can 
accelerate the flow to at most U rc f. Then, recalling £ ref = L rci /U rc {, we see that L ref = U 2 e{ /g. 
Recalling then that i? re f = p re f/(Pref g) and p re f = p rc f c^ rf , we see that L Te {/H Te f = 0(M 2 ). 
In this case the nondimensional momentum equation has the form 

+ V • (pUU) + ^ V(p - po) = -05 - p )^e r , 

which is consistent with the low Mach number nuclear burning model used in Bell et al. 
(2004). The assumption that t re f = U rc f/g is in fact unnecessarily restrictive; in a realistic 
physical scenario, even in the case of locally large density variations, the fluid accelerates with 
acceleration DU/Dt = a < g, and the relevant time scale would in fact be t re f = U rc {/a, i.e., 
the model is valid as the fluid accelerates, until the Mach number of the flow is no longer 
small. In other words, the assumption of small Mach number is sufficient to guarantee 
validity of the model. 

We note two important features of the model thus far. The first is that in both cases 
the perturbational pressure, which we will now denote as vr(x, t) = p(x, t) — po(r) satisfies 
n/po — 0(M 2 ). Thus in all but the momentum equation (where the 1/M 2 scaling requires 
the presence of Vn), we can substitute po for p. It is this approximation that decouples the 
pressure from the density in such a way as to filter acoustic waves from the solution. 

The second important feature is that the momentum equation can be retained in its 
original form, 

^ + V • (pUU) + Vvr = -Go - p )ge r , 

with no approximation to the buoyancy term and no assumption on the size of the pertur- 
bational density, as long as the actual acceleration of the flow is such that Mach number of 
the flow remains small. 

We now consider the implications of replacing p by po- The evolution of the non-reacting 
low Mach number system is described by the mass and momentum equations in combination 
with the enthalpy equation, but the system remains constrained by the equation of state 
(Eq. [5]), namely, p(p,X k ,T) = po(z). To complete the low Mach number model we re-pose 
the equation of state as a constraint on the velocity field, following closely the derivation in 
Bell et al. (2004) but retaining stratification effects. 

We begin by re-writing conservation of mass as an expression for the divergence of 
velocity: 
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Differentiating the equation of state (Eq. [5]) along particle paths, we can write 

DX k 



Dp dp 
15t = dp 



T,X k 



Dp dp 
~Dt + df 



or 



Dp 1 Dp 



p,x k 
DT 



E ' 



Dt 



p T 



Dt 



E 



dX, 



Dt 



Dt Pp y Dt 

with p p = dp/dp\ TXk , p T = dp/dT\ pXk , and p Xk = dp/dX k \ pT . 

We now require an expression for DT/Dt, which can be found by differentiating the 
enthalpy equation (Eq. [6]): 



Dh 



P- 



Dt 

or, gathering terms 



dh 



DT dh 

nX ~Dt + dp~ 

p,x k ^ 



T,X k 



Dt ^ 



dh 
dXl 



T,p 



DXk 
Dt 



TJt ■ 



DT 1 



Dt pc p 



Dp 
~D~t 



(9) 



where c p = dh/dT\ p Xk is the specific heat at constant pressure, = dh/dXk\ pT , and h p = 
dh/dp\ TXk for convenience. Substituting equation (9) into equation (8) and the resulting 
expression into equation (7) yields 



V- U 



= — f^(i-^)-i)^ + — f^-pE^*)+E^*) • 

PP P \PC P K P ' J Dt pp P \pc p ^ ^ k J 

Then, replacing p by Po{r), Dp/Dt becomes U ■ Vpo> an d, recalling the definition of H, the 
divergence constraint can be written 



V • U + all ■ Vp 



where we define 



Pt 
PPp \pc P 



( V ■ (kVT) - Pfa + Zh)uk ) + E 



a(p,T) = 



[1 - ph p )p T - pc p 
P 2 c P Pp 



Px k u k = S , (10) 



We note that for domains sufficiently smaller than a pressure scale height where Vpo can be 
neglected, equation (10) reduces exactly to the divergence constraint, equation (5) in Bell 
et al. (2004). 
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For the larger domains that are the target of this paper, we use the thermodynamic 
identities as outlined in Appendix A, to write 



1 

a 



where Fx = d(\ogp) / d(log p)\ s , and we have substituted po for p. 

In the case of terrestrial atmospheres and in the absence of compositional effects, I\ 
is replaced by the constant 7 = c p /c v , and using p = pRT with R the gas constant, this 
expression can be simplified to 

v.(p»= RH 



R/c p 



the pseudo-incompressible constraint as derived by Durran (1989). 

For stellar atmospheres the variation in Ti can be decomposed into two contributing 
factors: the background stratification and the local perturbation to that base state. For a 
nearly isentropically stratified base state and small perturbations to the base state, Ti is 

1 /T1 

close to constant, hence po oc p 10 . Neglecting expansion effects from thermal diffusion or 
reactions, this then reduces to the traditional anelastic approximation, 

V • (p U) = . 



For the more general case r lo = Ti(pQ,T ,p ) is not constant, but we can exploit the 
fact that both po and ri are functions only of r. It is straightforward (see Appendix B) to 
show that in this case the constraint can be written 

V ■ (J3o(r)U) = P S , (12) 

where 

/9b(r) = /3(0)exp( F (^) dr') . 

For the types of problems amenable to the low Mach number model, the density and 
temperature perturbations may be large, but even so the variation of Ti due to the pertur- 
bation is at most a few percent. For the examples in this paper we will neglect the local 
variation of Ti; this assumption will be re-examined in subsequent work. 

For the purposes of comparing the fundamental hydrodynamic behavior of the low Mach 
number model relative to established compressible and anelastic formulations, we will for the 
remainder of this paper neglect the effects of variation in composition, reactions, and thermal 
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conduction. Summarizing the low Mach number equation set for this specialized case and 
re-writing the momentum equation as an evolution equation for velocity instead, we now 
have 

dp 

dt 
dU 

~dt 

V ■ (foU) 

We note that this system contains three equations for three unknowns: density, velocity 
and pressure. The equation of state was used to derive the constraint thus to include it 
here would be redundant. When reactions and compositional effects are included in future 
work, the evolution equations for species and energy (in the form of temperature, entropy or 
enthalpy) will be added to this system, but for the simple hydrodynamical tests we present 
here this system is sufficient. 

3. Numerical Methodology for the Low Mach Number Model 

We discretize the low Mach number equation set derived in the previous section using an 
extension of the second-order accurate projection methodology developed for incompressible 
flows (Bell et al. 1989, 1991; Almgren et al. 1996; Bell & Marcus 1992; Almgren et al. 1998), 
and extended to low Mach number combustion (Pember et al. 1998; Day & Bell 2000) and 
to small-scale reacting flow for SNe la (Bell et al. 2004). We refer the reader to the above 
references for numerical examples demonstrating the second-order accuracy of the overall 
methodology, as well as for many of the details of projection methods. Here we present a 
brief overview of the numerical methodology as applied to this equation set. The absence 
of reactions and presence of flo in the projection steps are the key differences relative to the 
algorithm in Bell et al. (2004). 

In the projection approximation, explicit discretizations of the evolution equations are 
first used to approximate the velocity and thermodynamic variables at the new time, then an 
elliptic equation for pressure, derived from the constraint imposed on the new-time velocity, 
is solved to update the pressure and return the velocity field to the constraint surface. By 
contrast with the traditional anelastic approach, we have not replaced conservation of mass 
by the divergence constraint, which means we are able to evolve the density field with a 
conservative update, rather than invoking the equation of state to diagnose it. Thus the 
variables we update in each advection step are the velocity, density, and either temperature, 
enthalpy, or entropy. The low Mach number constraint constrains the evolution of the 



-V • (pU) , 

TT V7TT 1 V7 jP ~ gg) Q 

—U ■ \U V7T ge r 

P P 
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thermodynamic variables to the manifold defined by the equation of state. 

The discretization of the evolution equations is essentially a three-step process. First, 
we use an unsplit second-order Godunov procedure (Colella 1990) to predict a time-centered 
(t n+1 / 2 ) edge-based advection velocity, U ADY '*, using the cell-centered data at t n and the 
lagged pressure gradient from the interval centered at t n ~^ 2 . The provisional field, £/ ADV >* ; 
represents a normal velocity on cell edges analogous to a MAC-type staggered grid discretiza- 
tion of the Navier-Stokes equations (Harlow & Welch 1965). Figure 1 illustrates the MAC 
grid. However, f/ ADV >* fails to satisfy the time-centered divergence constraint (Eq. [12]). We 
apply a discrete projection by solving the elliptic equation 



J D MAC (^G MAC MAC ) = D MAC ((3 U ADY '*) - (3 S n+ K (13) 
p n 



for MAC , where £> MAC represents a centered approximation to a cell-based divergence from 
edge-based velocities, and G MAC represents a centered approximation to edge-based gradients 
from cell-centered data. The solution, </> MAC , is then used to define 

jjADY _ jjADV,* ^_^MAC iMAC 

jjKdv j g a seconc l_ orc l er accurate, staggered-grid vector field at t n+1 k that discretely satisfies 
the constraint (Eq. [12]), and is used for computing the time-explicit advective derivatives 
for U and p. 

We next explicitly update the density using a second-order accurate discretization of 
the mass equation. (We note here that this approach differs from both the anelastic equation 
set and the alternate form of the low Mach number equations as described in Lin 2005.) 

n+i _ -n \+\v7 ^rrAr>V\i n +y2 



p n+l = p " _ At [V ■ ( P U 

The final step of the integration procedure is to advance the velocity to the new time 
level. For this step we first obtain a provisional cell-centered velocity at t n+l using a time- 
lagged perturbational pressure gradient, 



P 



rm+1,* Tjn . 



+ p n+1 h [(U ADV ■ V)U] n+h = -Gn n ^ - (p"+ 1A - p )ge r 



At 

where p n+1 k = (p n + p n+1 )/2. At this point [/ n+1 >* does not satisfy the constraint. We 
apply an approximate projection to simultaneously update the perturbational pressure and 
to project U n+1 '* onto the constraint surface. In particular, we solve 

(TTn+l,* i \ a qn+l 
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for nodal values of <p, where is the standard bilinear finite element approximation to 
V- (A)/V)V with p evaluated at t n+ k. In this step, D is a discrete second-order operator that 
approximates the divergence at nodes from cell-centered data, and G = —D T approximates 
a cell-centered gradient from nodal data. (See Almgren et al. (1996) for a detailed discussion 
of this approximate projection; see Almgren et al. (2000) for a discussion of this particular 
form of the projection operand.) Finally, we determine the new-time cell-centered velocity 
field from 

and the new time-centered perturbational pressure from 

7T ?t+1/2 = . 

Specification of the initial value problem includes initial values for U and p at time t — 
and a description of the boundary conditions, but the perturbational pressure is not initially 
prescribed. To begin the calculation, then, the initial velocity field is first projected to ensure 
that it satisfies the divergence constraint at t — 0. Then initial iterations (typically two are 
sufficient) are performed to calculate an approximation to the perturbational pressure at 
t = At/2. 

In each step of the iteration we follow the procedure described above. In the first 
iteration we use 7r~ 1/2 = 0. At the end of each iteration we have calculated a new value of U 1 
and a pressure ir^ 2 . During the iteration procedure, we discard the value of U 1 , but define 
7T -1 / 2 = 7]- 1 / 2 . Once the iteration is completed, the above algorithm can be followed as written. 

4. Numerical Validation and Comparison 

4.1. Compressible Formulations 

We compare and contrast the low Mach number results with those obtained using two 
different discretizations of the fully compressible equation set, both implemented in the 
FLASH Code (Fryxell et al. 2000). The first is the piecewise parabolic method (PPM) (Colella 
& Woodward 1984), which is a high-order accurate, dimensionally split algorithm where the 
updates are done in one-dimensional sweeps, e.g. in two dimensions 

S^ 2 = X(At)Y(At)Y(At)X(At)S^ , 



where S = (p, pU, pE) is the state variable, X(At) is the operator that updates the state 
through At in time in the horizontal direction, and Y(At) updates the state by At in the 
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vertical direction. One PPM cycle updates the state through 2 At, switching the order of 
the directional operators midway through to retain second-order accuracy. PPM is the 
primary hydrodynamics algorithm used by the large-scale SNe la explosion modeling com- 
munity (Ropke & Hillebrandt 2005; Plewa et al. 2004; Gamezo et al. 2003). The FLASH 
implementation of PPM has been well validated (Calder et al. 2002), and serves as a good 
basis for comparisons with the low Mach number algorithm. 

A numerical issue that arises in fully compressible simulations, but not with the low 
Mach number approach, is the difficulty of maintaining a quiet hydrostatic atmosphere. 
Small displacements from hydrostatic equilibrium (HSE) can generate sound waves through- 
out the atmosphere, which, if unchecked, can lead to ambient velocities that can swamp 
the process being studied. The hydrostatic equilibrium improvements described in Zingale 
et al. (2002), which remove the hydrostatic pressure from the pressure jump across the inter- 
faces in the Riemann problem, were used for all PPM runs. The upper and lower boundary 
conditions are the hydrostatic boundaries described in that same paper, with the pressure 
and density modified according to hydrostatic equilibrium, and the velocities given a zero 
gradient. 

The second compressible algorithm we consider is a second-order unsplit method fol- 
lowing Colella (1990). At low Mach number, dimensionally split methods can have trouble 
producing realistic velocity fields, as will be shown in the bubble rise comparison. In the 
unsplit formulation, the cell averages are updated in all directions at once. A critical part 
of the unsplit method is that the interface reconstructions contain a transverse flux term 
that explicitly couples in the information from the corner cells. Fourth-order accurate slope 
limiting is used in the central differences, as well as high-order reconstruction of the states for 
the transverse Riemann problem, as described in (Colella 1985). This is the same procedure 
used in predicting the interface states in the low Mach number method presented here. This 
method was extended to handle general equations of state following the procedure in Colella 
& Glaz (1985), adding an additional transverse flux piece to the interface reconstruction of 
7, to be consistent with the unsplit reconstructions. This was put into the FLASH frame- 
work for the present simulations. Both the PPM and unsplit solvers use the same two-shock 
Riemann solver described in Colella & Glaz (1985). For both the split and unsplit solvers, a 
CFL number of 0.8 was used based on the sound speed. 

4.2. Anelastic Approach 

The low Mach number equations and the anelastic equation set are derived differently. 
Both equation sets assume a low Mach number, equivalently a small pressure perturbation 
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from the background state. However, the anelastic equation set assumes both small den- 
sity and small temperature perturbations as well. As noted earlier the velocity constraints 
resulting from these two derivations are strikingly similar, and in fact equivalent for an isen- 
tropically stratified background state. Even in the non-isentropic background considered 
here, the differences between p and (3q are small. 

However, because the anelastic approximation assumes small density and temperature 
perturbations, approximations are made to the buoyancy term in the momentum equation. 
These follow from the observation that since the perturbational density was neglected in 
the continuity equation in order to derive the velocity constraint, the continuity equation 
cannot be used to evolve the perturbational density. Thus an alternative formulation of the 
buoyancy term must be used. A typical anelastic model evolves temperature or entropy, 
and constructs the buoyancy forcing term from that field using a linearized approximation. 
Following the derivation by Braginsky & Roberts (1995) that combines parts of the pressure 
gradient and buoyancy terms, we consider the following form of the anelastic equations: 

DU (p'\ fdpo\ „, 

^ = », 

Dt 

V • (poU) = . 

Here S is entropy, S' = S—So, p' = p—po, and we have neglected viscosity, thermal diffusivity 
and the gravitational potential perturbation. We note that in the case of small density 
and temperature perturbations, simulations using the low Mach number equations and the 
anelastic approximation give indistinguishable results, thus for the numerical comparisons we 
focus on problems with finite amplitude perturbations as described in the next subsection. 

4.3. Bubble Rise Comparison 

We present three sets of two-dimensional calculations of a rising bubble in a stellar 
environment. The one- dimensional background state (po, T , p ) is calculated using the 
Kepler code (Weaver et al. 1978), to evolve a Chandrasekhar mass white dwarf until the 
central temperature reaches 7 x 10 8 K, representing conditions just before ignition. We map 
a portion of the one-dimensional model onto a uniform two-dimensional grid, and place it into 
hydrostatic equilibrium with a constant gravitational acceleration (g = — 1.9 x 10 10 cm s~ 2 ). 
We further simplify by ignoring metric terms associated with the radial coordinate and view 
the domain as Cartesian. We note that neither the constant gravity assumption nor the 
simplified metric is a limitation of any of the methods presented here, but is chosen in these 
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comparison simulations for simplicity. The density structure of the model is illustrated in 
Figure 2. 

All bubbles begin in pressure equilibrium with the background state and are defined by 
a simple temperature perturbation, from which the density perturbation is calculated. We 
consider three different cases, which we will distinguish by the maximum temperature at the 
center of the bubble, T max . The temperature profile of the bubble is then defined by 



and (x cent , r cent ) = (2.5 x 10 7 , 6.25 x 10 7 ) cm, 5 = 1.25 x 10 6 cm in a domain from x = cm to 
5 x 10 7 cm and r = 5 x 10 cm to 10 8 cm. The stellar equation of state is then used to compute 
p given T and po- This profile was chosen to give a smooth transition from the ambient 
temperature to the perturbed temperature, thus minimizing the effects of the numerical 
slope limiters present in the different hydrodynamics methods. Due to the short time scale 
of the problem thermal diffusivity is neglected. For all bubble calculations presented the 
grid has uniform resolution of 384 x 384; the adaptive gridding features of all the codes are 
turned off. 

Figures 3 and 6 present comparisons of simulations using the low Mach number ap- 
proach, two different discretizations of the fully compressible equation set, as well as the 
anelastic and incompressible equation sets. The low Mach number, anelastic and incom- 
pressible results are calculated using the projection method approach described in the pre- 
vious section. The only differences in methodology occur in the coefficient of velocity in the 
projection, and in the construction of the buoyancy term. Each of these methods was run 
at a CFL number of 0.9 based on the maximum advection velocity. 

Figure 3 shows the temperature evolution for T max = 6 x 10 9 K. This corresponds to 
an Atwood number for the bubble of approximately 0.079. In this simulation, the bubble 
reaches a Mach number of about 0.2. In addition to the PPM and unsplit compressible 
solvers, traditional anelastic and incompressible solvers are shown for comparison. The 
low Mach number method closely tracks the two compressible solvers. The incompressible 
and anelastic results demonstrate the effects of their respective assumptions. The velocity 
constraint for the anelastic model is sufficiently accurate to capture the bubble rise, but 
because of the linearity of the density-temperature relationship in the anelastic approach, 
the buoyancy term is too small in the anelastic simulation. By contrast, the incompressible 
simulation contains the full buoyancy term but due to the incompressibility constraint the 




where 
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bubble cannot expand; consequently, it reaches neutral buoyancy at a much lower level and 
stops rising. We do not follow the bubbles past the point where nonlinear instabilities along 
the sides begin to dominate the evolution. 

A more detailed comparison of the results from Figure 3 is provided in Figure 4, where 
temperature contours of the low Mach number solution are superimposed on temperature 
contours of the the unsplit and PPM solutions. Here we see a large degree of overlap, 
demonstrating that the bubbles have the same rise velocity and size independent of the 
algorithm. Figure 5 shows the Mach number of the PPM and low Mach number methods, 
further demonstrating the agreement between the two sets of results, with the exception of 
the unphysical loss of symmetry in the PPM simulation. 

Figure 6 shows the temperature evolution for the lower peak temperature case, T max = 
1 x 10 9 K, and corresponding Atwood number of 0.0024. Over the course of this comparison, 
the Mach number remains below 0.05. Again, we observe the agreement between the low 
Mach number and compressible solvers, with the exception of the late-time breakdown of 
the PPM solution, indicated by the large amplitude temperature oscillations dominating the 
flow behind the bubble. These oscillations reflect the poor performance of operator split 
algorithms for very low speed flows. 

A more detailed comparison of the results from Figure 6 is shown in Figure 7, again 
superimposing temperature contours from the low Mach number and compressible formula- 
tions. We again see good agreement, with the exception of the breakdown of PPM at late 
times. 

Timings of the PPM and low Mach number code were made on a single processor 
(1.53 GHz Athlon MP) using the Intel 8.1 compilers. Both codes were compiled with the 
same compiler optimization flags: -03 -ip -ipo. FLASH was set to run with 16 x 16 zone 
blocks, instead of the default 8x8 for better performance with uniform gridding. The 
6 x 10 9 K bubble required 2148 timesteps, taking 14200 s to evolve the bubble to 0.25 s in 
simulation time. By comparison, the low Mach number solver took 246 timesteps and 1480 s, 
about an order of magnitude speed-up. For the 10 9 K bubble, the PPM solver took 7842 
timesteps, taking 52100 s to evolve the bubble for 1 s of simulation time, while the low Mach 
number solver took 252 timesteps and 1560 s. As expected, the performance gap increases as 
the Mach number decreases. The unsplit compressible algorithm takes approximately twice 
as much time to run as PPM, primarily due to the additional transverse Riemann solves 
required. 

Finally, in Figures 8 and 9 we compare the low Mach number and anelastic models for a 
bubble with T max = 3.5 x 10 8 K. This regime is inaccessible to the compressible formulation, 
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with a peak Mach number during the calculation of 0.012. We note that, as expected, as 
the Atwood number decreases the fidelity of the anelastic approximation improves. 

The three cases presented in this section demonstrate successful application of the low 
Mach number approach, as well as failure of the compressible approach for low-speed flows, 
and failure of the anelastic approach for flows with large density or temperature variations. 
The low Mach number approach, like the other methods, has limits to its applicability. 
Specifically, as the flow speed, hence the Mach number, increases, typically the thermody- 
namic pressure will diverge from po and the assumptions underlying the low Mach number 
approach will be violated. As with the anelastic model, numerical simulations using the 
low Mach number model will continue to yield what appear to be reasonable solutions even 
as the underlying assumptions are violated, but the solutions will no longer be physically 
relevant. 

In the case of a bubble rise without heat sources, it is difficult to numerically demonstrate 
the failure of the low Mach number method even as the Mach number increases. However, 
the divergence of the low Mach number model from the compressible solution in the case of 
large Mach number will be discussed more thoroughly in the subsequent paper that discusses 
the behavior of the low Mach number model in the presence of heating. 

5. Conclusions 

We have introduced a new method for following low speed, stratified flows in astrophys- 
ical conditions and have demonstrated, through comparison with compressible and anelastic 
codes, that this algorithm performs well in the range of Mach number from near zero to about 
0.2. The increased computational efficiency associated with a low Mach number formulation 
makes it an ideal tool for investigations of the convective/ignition phase of SNe la. However, 
to be applicable in this setting, a number of generalizations to the methodology will need to 
be developed. In particular, we will need to extend the the algorithm to include the effects 
of variation in composition, reactions and thermal conduction. In addition, once a flame is 
established it will be necessary to include sub-grid models for turbulent flame propagation 
that will enable the methodology to be used to simulate the evolution of the early phases of 
the explosion. These issues will be addressed in future work. 

We thank Gary Glatzmaier, Michael Kuhlen, and Tami Rogers for helpful discussions 
regarding the anelastic hydrodynamics method, and Alan Calder and Jonathan Dursi for 
helpful comments on the manuscript. We especially thank Stan Woosley for numerous discus- 
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DOE grant No. DE-FC02-01ER41176 to the Supernova Science Center/UCSC, and by the 
NASA Theory Program (NAGW-5- 12036). 



A. Simplification of a 



In this appendix we derive a simplified expression for a introduced in equation 11. 
We refer to Chapter 9 of Cox & Giuli (1968) (CG, in this appendix), for a number of 
thermodynamic identities. 



We begin by rewriting h p , using 

dh 
dp 
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h p Pp 
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Therefore, 
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Putting this into a, we have 
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For a generalized equation of state, there are three principal adiabatic exponents which 
relate the various differentials (dp, dT, and dp). For an ideal gas, they are all equivalent. 
Here, we use Ti (CG eq. 9.88): 

Y _ f d^-P^ 



This is related to the ratio of specific heats, 7 via 
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(CG eq. 9.87) where 



dlnp 
dlnp 



P. 
V 



Pp 



(CG eq. 9.82) is the "density exponent in the pressure equation of state." For an ideal gas, 
X P = 1, and Ti = 7. Taking all of this together, we see that 



PPp Tip 



Putting this into our a expression, 
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Motivated by the ideal gas result that a = 1/(72?), we want to show that the quantity in the 
square brackets in equation (Al) reduces to cy. 



The specific heats are related by (CG eq. 9.84) 

E ( dlnE\ xt P Xt 
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T V dlnp / T x P pT x P 
The temperature exponent is defined as 
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(CG eq. 9.81), so equation (A2) simplifies to 
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which substitutes directly into equation (Al) to yield 
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We note that Ti varies slowly throughout the white dwarf and is a quantity that is already 
returned by the tabular equation of state. 
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B. Derivation of (3 



We seek a function (3(z) such that 



L -V ■ (PU) = (V ■ U) + -*-[/ ■ Vpo 



(3(z) Tip 

We expand V ■ (/3U) — (3{SJ ■ U) + U • V ' (3 and note that for the equality to hold we would 
need 

U-Vf3 = —U ■ Vp , 



/3(z) Tipo 

1 1 

- w (3' = - — wp 

Since we want this to hold for all w, we are left with 

i^ = A_ . 

P 

We integrate this up from z = : 

f ^ , „ - , >r 
Jo P Jo dz Jo r iPo 



so 

n y 

Po 



ln([3(z)) - ln(/3(0)) = / (-^) ^ 

7o r iPo 

or 

0(z) = 0(0) exp ( f\^)dz' 
\Jo r iPo 

We note that this also can be written as the recursive relationship 

P(z k ) = P(Zk-i) exp I / {--^-)dz' 



.! F lP0' 



exploiting the hydrostatic equilibrium of the base state. This equation is the one we use to 
numerically compute P(z); we let (3(0) = po(0). 
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Fig. 1. — Illustration of the MAC-type grid showing the advective velocities (w ADV , v ADV ) 
and (j) MAC for the cell. An unsplit Godunov method is used to predict these advective 
velocities from the cell centered velocities (wjj, Vij) in the surrounding cells. 
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Fig. 2. — Initial model generated by the Kepler code (dotted) and uniformly gridded, con- 
stant gravity portion used in the bubble simulations (solid). 
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Fig. 3. — Bubble evolution for five different algorithms. Here, the peak temperature is 
6 x 10 9 K. We see good agreement between the two compressible codes (PPM and unsplit) 
and the low Mach number algorithm. 
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Fig. 4. — Detailed comparison of the temperature field for PPM (blue), the unsplit com- 
pressible algorithm (red), and the low Mach number algorithm (green), for the 6 x 10 9 K 
perturbation. Contours are shown at 3 x 10 9 and 4.5 x 10 9 K. 
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Fig. 5. — Mach number comparison for the 6 x 10 9 K bubble. At early times (a.), the PPM 
results show a sound wave from the initial perturbation just about to exit through the top 
of the domain. This is not present in the low Mach number it filters sound waves. 

At late times (b.), the flow has attained a Mach number of almost 0.2 in places. 
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Fig. 6. — Bubble evolution for four different algorithms. Here, the peak temperature is 
1 x 10 9 K. As with the 6 x 10 9 K case, the PPM, unsplit, and low Mach number results are 
in good agreement. 
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Fig. 7. — Detailed comparison of the temperature field for PPM (blue), the unsplit com- 
pressible algorithm (red), and the low Mach number algorithm (green), for the 1 x 10 9 K 
perturbation. Contours are shown at 5 x 10 8 and 7.5 x 10 8 K. 




Fig. 8. — Bubble evolution the low Mach number and anelastic algorithms, showing good 
agreement for the two different methods.. Here, the peak temperature is 3.5 x fO 8 K. 
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Fig. 9. — Detailed comparison of the temperature field for the anelastic algorithm (purple) 
and the low Mach number algorithm (green), for the 3.5 x 10 8 K perturbation. The contours 
are at 2.5 x 10 8 and 3 x 10 8 K 



